clear al; clc; close all;
save_res=1;

alpha_s=deg2rad(30);%1
beta_s=deg2rad(0);

beta_j=deg2rad(0);

alpha = deg2rad(-180:2:180);
beta = deg2rad(-90:3:90);

Fa = ones(length(beta),1)*(1+cos(alpha-pi/2)).^2;
[alpha_mesh,beta_mesh]=meshgrid(alpha,beta);
[xa,ya,za]=sph2cart(beta_mesh,alpha_mesh,Fa);

figure(1);
surf(xa,ya,za);
xlabel('X');ylabel('Y');zlabel('Z');
title('Radiation pattern for one AM');
axis equal

F=nan(length(beta),length(alpha));

q_dB = 3; %1
q = 10^(q_dB/10);

figure(2)
pos = get(gcf, 'Position'); pos(3) = 800; set(gcf, 'Position',pos);

for alpha_j=deg2rad(30:5:95)
    
    C=H(alpha_j, beta_j);
    Dn=q*C*C'+eye(7);
    Hf=H(alpha_s,beta_s);
    beta_w=Dn\Hf/(Hf'*(Dn\Hf));
    
    for a=1:length(alpha)
        for b=1:length(beta)
            U=beta_w'*H(alpha(a),beta(b));
            F(b,a)=abs(U)^2;
        end
    end
   F=Fa.*F;
   [x,y,z]=sph2cart(beta_mesh,alpha_mesh,F);
   
   b0=ceil(length(beta)/2);
   Fb0=F(b0,:);
   
   subplot(1,2,1)
   surf(x,y,z)
   xlabel('X');ylabel('Y');zlabel('Z');
   axis equal
   
   subplot(1,2,2)
   polar(alpha,Fb0);
   hold on
   polar([alpha_s alpha_s],[0 max(Fb0)],'g');
   polar([alpha_j alpha_j],[0 max(Fb0)],'r');
   hold off
   drawnow
   
   if save_res
       s=sprintf('pic/DN_alpha_j_%03.0f.png',round(rad2deg(alpha_j)));
       saveas(gcf,s,'png');
       fprintf('Figure is saved at %s\n',s)
   end
end